/*
 * Created on 26-lug-2005
 *
 * To change the template for this generated file go to
 * Window&gt;Preferences&gt;Java&gt;Code Generation&gt;Code and Comments
 */
package it.unibo.deis.lia.filters;

import java.util.Vector;

public class Fourier implements IFilter{
	//private boolean verbose=true;
	
	public IFilter genFilter(){
		return new Fourier();
	}
	
	public Double nextValue(Vector<Double> valueSet){//, double costante) {
		// pmax in the maximum fourier coefficient index
		// exploited in the inverse transform
		int pmax=1; // pmax < valori.size()/2 (io penso)
		if(valueSet.size()==0)return null;
		if(valueSet.size()/2<pmax)return valueSet.lastElement();
		
		/*%
		 % DFT in a direct implementation
		% Enter Data in y
		y=[7.6 7.4 8.2 9.2 10.2 11.5 12.4 13.4 13.7 11.8 10.1 ...
		    9.0 8.9 9.5 10.6 11.4 12.9 12.7 13.9 14.2 13.5 11.4 10.9 8.1];

		% Get length of data vector or number of samples
		N=length(y);*/
		
		
		/*% Compute Fourier Coefficients
		for p=1:N/2+1
		    A(p)=0;
		    B(p)=0;
		    for n=1:N
		        A(p)=A(p)+2/N*y(n)*cos(2*pi*(p-1)*n/N)';
		        B(p)=B(p)+2/N*y(n)*sin(2*pi*(p-1)*n/N)';
		    end   
		end

		A(N/2+1)=A(N/2+1)/2;*/
		
		Vector<Double> A=new Vector<Double>();
		Vector<Double> B=new Vector<Double>();
		int N=valueSet.size();
		//if(verbose)System.out.println("N "+N);
		
		//double maxMod=0;
		for(int p=0;p<N/2+1;p++){
			double tempA=0;
			double tempB=0;
			for(int n=0;n<N;n++){
				//System.out.print(n+"\t");
				tempA=tempA+2.0/N*valueSet.elementAt(n)*Math.cos(2.0*Math.PI*p*n/N);
				tempB=tempB+2.0/N*valueSet.elementAt(n)*Math.sin(2.0*Math.PI*p*n/N);
			}
			//System.out.println();
			
			/*if(verbose){//&& p==1){
				//System.out.println("Armonic:"+p+"        tempA:"+tempA+" tempB:"+tempB);
				double real=tempA/2;
				double imm=-tempB/2;
				double mod=Math.sqrt(real*real+imm*imm);
				double ph=Math.atan(imm/real);
				
				//if(p>=1&&p<=3)modCount=modCount+mod;
				if(p>=1){//&&p<=3){
					if(maxMod<mod)maxMod=mod;
				}
				//System.out.println("Armonic:"+p+"        mod:"+mod+" ph:"+ph+" mod*ph:"+mod*ph);
			}*/
			
			A.addElement(tempA);
			B.addElement(tempB);
		}
		//System.out.println("maxMod: "+maxMod);
		//A(N/2+1)=A(N/2+1)/2; ???
		A.setElementAt((A.elementAt(N/2)).doubleValue()/2,N/2);
		//if(verbose)System.out.println("A "+A);
		//if(verbose)System.out.println("B "+B);

		/*% Reconstruct Signal - pmax is number of frequencies used in increasing order
		pmax=13;

		for n=1:N
		    ynew(n)=A(1)/2;
		    for p=2:pmax
		        ynew(n)=ynew(n)+A(p)*cos(2*pi*(p-1)*n/N)+B(p)*sin(2*pi*(p-1)*n/N);
		    end
		end   */

		Vector<Double> ynew=new Vector<Double>();
		for(int n=0;n<N;n++){
			double ynewTemp=A.elementAt(0)/2;
			for(int p=1;p<pmax;p++){
				//if(verbose)System.out.println("p "+p);
				//if(verbose)System.out.println("A.size() "+A.size());
				ynewTemp=ynewTemp
					+A.elementAt(p)*Math.cos(2*Math.PI*p*n/N)
					+B.elementAt(p)*Math.sin(2*Math.PI*p*n/N);
			}
			ynew.addElement(ynewTemp);
			//if(verbose)System.out.println("ynewTemp "+ynewTemp);
		}
		
		return ynew.lastElement();

	}
}
		
		
		
		
		
/*----------------------------------------------------------------------------------
		Compact  Implementation
		This implementation uses matrix products and is described in Appendix H.

		% DFT in a compact implementation

		% Enter Data in y
		y=[7.6 7.4 8.2 9.2 10.2 11.5 12.4 13.4 13.7 11.8 10.1 ...
		    9.0 8.9 9.5 10.6 11.4 12.9 12.7 13.9 14.2 13.5 11.4 10.9 8.1];

		N=length(y);		 

		% Compute the matrices of trigonometric functions
		p=1:N/2+1;
		n=1:N;

		C=cos(2*pi*n'*(p-1)/N);
		S=sin(2*pi*n'*(p-1)/N);

		% Compute Fourier Coefficients
		A=2/N*y*C;
		B=2/N*y*S
		A(N/2+1)=A(N/2+1)/2;		 

		% Reconstruct Signal - pmax is number of frequencies used in increasing order
		pmax=13;
		ynew=A(1)/2+C(:,2:pmax)*A(2:pmax)'+S(:,2:pmax)*B(2:pmax)';

		*/

		/*double realarr[] = new double[100];
		double imagarr[] = new double[100];
		double diff[] = new double[100];
/*----------------------------------------------------------------------------------
	 	
		// transform
		int i=valueSet.size();
		for (int j = 0; j<i; j++)                      //loop for frequency index
		{
		  double real = 0.0;
		  double imag = 0.0;                       //clear variables
		  for(int k = 0; k<i; k++)                     //loop for sums
		  {
		    real += valueSet.elementAt(k)*Math.cos((2*Math.PI*k*j)/i);
		    imag += valueSet.elementAt(k)*Math.sin((2*Math.PI*k*j)/i);
		  }
		 
		  Integer total = new Integer(i);          //need to get doublevalue of i
		  double totald = total.doubleValue();
		  realarr[j] = real/totald;
		  imagarr[j] = imag/totald;
		  
		  //s = Format.sprintf( "%d\t%f\t%f\n", p.add(j).add(real/totald).add(imag/totald));
//		  System.out.println("*** "+j+" "+real/totald+" "+imag/totald);
		  //s = Format.sprintf( "%d\t%f\n", p.add(j).add((Math.sqrt(real*real + imag*imag))/totald) );
		  System.out.println("*** "+j+" "+(Math.sqrt(real*real + imag*imag))/totald);
		  //out.write(s);
		}
		//out.close();}/**/
		
		/* //inverse transform
		
		for (int j=0; j<i; j++)
		{
		  double real1 = 0.0
		  ; double imag1 = 0.0;
		  for(int k=0; k<i; k++)
		  {
		    real1 += (realarr[k] * Math.cos(2*Math.PI*k*j/i)) + (imagarr[k] * Math.sin(2*Math.PI*k*j/i));
		    imag1 += (imagarr[k] * Math.cos(2*Math.PI*k*j/i)) - (realarr[k] * Math.sin(2*Math.PI*k*j/i));
		  }
		  
		  diff[j] = Math.abs(valueSet.elementAt(j) - real1);
		  //s1 = Format.sprintf( "%d\t%e\n", p.add(j).add(diff[j]) );
//		  System.out.println("*** "+j+" "+diff[j]);
		  //diffout.write(s1);
		  //s = Format.sprintf( "%d\t%f\t%f\n", p.add(j).add(real1).add(imag1) );
//		  System.out.println("*** "+j+" "+real1+" "+imag1);
		  //inverseout.write(s);
		  
		}/**/
		
